{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 实验八 矩阵特征值及特征向量计算\n",
    "\n",
    "- 正幂法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def power_iteration(A, m0=1, u0=None, eps=1e-8, max_steps=500, verbose=False):\n",
    "    \"\"\"正幂法（power iteration a.k.a. the power method）计算矩阵 A 的按模最大特征值、特征向量\n",
    "    \n",
    "    Args:\n",
    "        A:   np_array_like 待求特征值的矩阵 (nxn)\n",
    "        m0:  float 初始特征值\n",
    "                default m0=1\n",
    "        u0: np_array_like 初始特征向量（n）\n",
    "                default u0=None: 取 u0 = (1, 1, ..., 1)\n",
    "        eps: float 精度要求\n",
    "                default eps=1e-8\n",
    "        max_steps: int 最大迭代次数\n",
    "                default max_steps=1000\n",
    "        verbose: bool, 若为 True 则打印出每一步的结果\n",
    "                default verbose=False\n",
    "        \n",
    "    Returns:\n",
    "        (m, u, k): 在 max_steps 次迭代以内得到第一组的满足 eps 的结果\n",
    "            m: float 所求主特征值\n",
    "            u: np.array 相应的特征向量\n",
    "            k: int 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 u0 尺寸不匹配\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "        \n",
    "    _shape = np.shape(A)\n",
    "    if len(_shape) < 2 or _shape[0] != _shape[1]:\n",
    "        raise ValueError(f\"unexpected A, shape: {_shape}\")\n",
    "    \n",
    "    if u0 is not None:\n",
    "        if len(u0) != _shape[0]:\n",
    "            raise ValueError(f\"A (shape={_shape}) and u0 (len={len(u0)}) not match\")\n",
    "    else:  # not u0\n",
    "        u0 = np.ones(_shape[0])\n",
    "        \n",
    "    m = m0\n",
    "    u = u0\n",
    "\n",
    "    for k in range(int(max_steps)):\n",
    "        if verbose:\n",
    "            print(k, m, u)\n",
    "        \n",
    "        m_prev = m\n",
    "        \n",
    "        v = np.dot(A, u)\n",
    "        mi = np.argmax(np.abs(v))\n",
    "        m = v[mi]\n",
    "        u = v / m\n",
    "        \n",
    "        if abs(m - m_prev) <= eps:\n",
    "            break\n",
    "            \n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). \"\n",
    "                        f\"The last result: m = {m}, u={u}\")\n",
    "        \n",
    "    if verbose:\n",
    "        print('result of power_iteration:', m, u, k+1)\n",
    "\n",
    "    return m, u, k+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "算一个课本上的例子（P159 例1）："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1 [1. 1. 1.]\n",
      "1 17.0 [0.41176471 1.         0.58823529]\n",
      "2 9.470588235294118 [0.52795031 1.         0.82608696]\n",
      "3 11.58385093167702 [0.49276139 1.         0.72600536]\n",
      "4 10.831635388739945 [0.50200485 1.         0.75743775]\n",
      "5 11.049799514875502 [0.49945569 1.         0.74783731]\n",
      "6 10.985906091381928 [0.50014864 1.         0.75061668]\n",
      "7 11.003953118800313 [0.49995948 1.         0.74982713]\n",
      "8 10.998903290037243 [0.50001105 1.         0.75004801]\n",
      "9 11.000302582696586 [0.49999699 1.         0.74998675]\n",
      "10 10.999916852432536 [0.50000082 1.         0.75000364]\n",
      "11 11.000022790833176 [0.49999978 1.         0.749999  ]\n",
      "12 10.999993763594336 [0.50000006 1.         0.75000027]\n",
      "13 11.000001704609195 [0.49999998 1.         0.74999993]\n",
      "14 10.999999534421143 [0.5        1.         0.75000002]\n",
      "15 11.000000127100696 [0.5        1.         0.74999999]\n",
      "16 10.999999965313513 [0.5  1.   0.75]\n",
      "17 11.00000000946407 [0.5  1.   0.75]\n",
      "18 10.999999997418142 [0.5  1.   0.75]\n",
      "result of power_iteration: 11.000000000704278 [0.5  1.   0.75] 19\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(11.000000000704278, array([0.5 , 1.  , 0.75]), 19)"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "power_iteration(A, eps=1e-8, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "用随机矩阵测试，和 numpy 提供的特征值算法比较："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "by power_iteration: 50.445839029418046\n",
      "by   np.linalg.eig: (50.44583922496655+0j)\n"
     ]
    }
   ],
   "source": [
    "A = np.random.random((100, 100))\n",
    "# print(A)\n",
    "\n",
    "r = power_iteration(A, eps=1e-5)\n",
    "print('by power_iteration:', r[0])\n",
    "\n",
    "a = np.linalg.eig(A)[0]\n",
    "ar = [x for x in a if not np.iscomplex(x)]\n",
    "i = np.argmax(np.abs(ar))\n",
    "print('by   np.linalg.eig:', ar[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 反幂法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse_iteration(A, m0=1, u0=None, eps=1e-8, max_steps=500, verbose=False):\n",
    "    \"\"\"反幂法（inverse iteration a.k.a. inverse power method）计算矩阵 A 的按模最小特征值、特征向量\n",
    "    \n",
    "    Args:\n",
    "        A:   np_array_like 待求特征值的矩阵 (nxn)\n",
    "        m0:  float 初始特征值\n",
    "                default m0=1\n",
    "        u0: np_array_like 初始特征向量（n）：要求无穷范式=1，通常取 u0 = (1, 1, ..., 1)\n",
    "                default u0=None: 取 u0 = (1, 1, ..., 1)\n",
    "        eps: float 精度要求\n",
    "                default eps=1e-8\n",
    "        max_steps: int 最大迭代次数\n",
    "                default max_steps=1000\n",
    "        verbose: bool, 若为 True 则打印出每一步的结果\n",
    "                default verbose=False\n",
    "        \n",
    "    Returns:\n",
    "        (m, u, k): 在 max_steps 次迭代以内得到第一组的满足 eps 的结果\n",
    "            m: float 所求主特征值\n",
    "            u: np.array 相应的特征向量\n",
    "            k: int 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 u0 尺寸不匹配\n",
    "                    或 u0 = 0\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "        \n",
    "    _shape = np.shape(A)\n",
    "    if len(_shape) < 2 or _shape[0] != _shape[1]:\n",
    "        raise ValueError(f\"unexpected A, shape: {_shape}\")\n",
    "    \n",
    "    if u0 is not None:\n",
    "        if len(u0) != _shape[0]:\n",
    "            raise ValueError(f\"A (shape={_shape}) and u0 (len={len(u0)}) not match\")\n",
    "        if np.all(u0 == 0):\n",
    "            raise ValueError(f'bad u0: u0 == 0')\n",
    "    else:  # not u0\n",
    "        u0 = np.ones(_shape[0])\n",
    "        \n",
    "    m = m0\n",
    "    u = u0\n",
    "\n",
    "    for k in range(int(max_steps)):\n",
    "        if verbose:\n",
    "            print(k, 1/m, u)\n",
    "        \n",
    "        m_prev = m\n",
    "        \n",
    "        v = np.linalg.solve(A, u)\n",
    "        mi = np.argmax(np.abs(v))\n",
    "        m = v[mi]\n",
    "        u = v / m\n",
    "                \n",
    "        if abs(m - m_prev) <= eps:\n",
    "            break\n",
    "        \n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). \"\n",
    "                        f\"The last result: 1/m = {1/m}, u={u}\")\n",
    "    \n",
    "    if verbose:\n",
    "        print('result of inverse_iteration:', 1/m, u, k+1)\n",
    "        \n",
    "    return 1/m, u, k+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1.0 [1. 1. 1.]\n",
      "1 2.75 [-0.25        0.41666667  1.        ]\n",
      "2 -1.7368421052631582 [-0.39473684 -0.25877193  1.        ]\n",
      "3 -1.522306525037936 [-0.27587253 -0.28244815  1.        ]\n",
      "4 -1.773137541073161 [-0.24854257 -0.33791831  1.        ]\n",
      "5 -1.850713802809302 [-0.22932389 -0.36045702  1.        ]\n",
      "6 -1.9067348353996336 [-0.21874548 -0.37508306  1.        ]\n",
      "7 -1.9394212652843035 [-0.21209945 -0.38385382  1.        ]\n",
      "8 -1.9604633754301686 [-0.20791023 -0.38945545  1.        ]\n",
      "9 -1.973976092920489 [-0.20520426 -0.39306055  1.        ]\n",
      "10 -1.9828015632828082 [-0.20343978 -0.3954137   1.        ]\n",
      "11 -1.988599435522524 [-0.2022801  -0.39695986  1.        ]\n",
      "12 -1.9924284508421255 [-0.20151431 -0.39798092  1.        ]\n",
      "13 -1.9949649983827589 [-0.201007   -0.39865733  1.        ]\n",
      "14 -1.996648958194817 [-0.20067021 -0.39910639  1.        ]\n",
      "15 -1.997768464464223 [-0.20044631 -0.39940492  1.        ]\n",
      "16 -1.998513415490028 [-0.20029732 -0.39960358  1.        ]\n",
      "17 -1.999009434502497 [-0.20019811 -0.39973585  1.        ]\n",
      "18 -1.9993398409805123 [-0.20013203 -0.39982396  1.        ]\n",
      "19 -1.9995599908120016 [-0.200088   -0.39988266  1.        ]\n",
      "20 -1.9997067035591078 [-0.20005866 -0.39992179  1.        ]\n",
      "21 -1.9998044881537032 [-0.2000391  -0.39994786  1.        ]\n",
      "22 -1.999869667263002 [-0.20002607 -0.39996524  1.        ]\n",
      "23 -1.9999131152833085 [-0.20001738 -0.39997683  1.        ]\n",
      "result of inverse_iteration: -1.999942078533036 [-0.20001158 -0.39998455  1.        ] 24\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-1.999942078533036, array([-0.20001158, -0.39998455,  1.        ]), 24)"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "inverse_iteration(A, eps=1e-5, max_steps=40, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "基于实验 6 写的 LU 分解:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lu 分解函数，来自 ex6/src/lu.py\n",
    "\n",
    "def lu(a, sequence=False, swap_times: list = {}):\n",
    "    \"\"\"\n",
    "    「高斯消去法」的 LU 分解.\n",
    "\n",
    "    本函数可以计算「列主元高斯消元法」、「顺序高斯消元法」的 LU 分解，\n",
    "    通过参数 sequence 控制，默认 sequence=False 使用「列主元高斯消元法」。\n",
    "\n",
    "    Args:\n",
    "        a: np_array_like 方阵 (nxn)\n",
    "        sequence: bool, True 则使用顺序高斯消去法，False 为列主元的高斯消去法\n",
    "            default: sequence=False\n",
    "        swap_times: 这是一个**输出**用的变量，只有传入 dict 变量时才有效。\n",
    "            若使用「列主元高斯消元法」（sequence=False）\n",
    "            则，置 swap_times['swap_times'] = 行交换次数。\n",
    "            这个值正常的输出中不需要，但在一些问题，比如，\n",
    "            利用 LU 分解求行列式时，得到 swap_times 会很有帮助。\n",
    "\n",
    "    Returns:\n",
    "        (l, u, p): result\n",
    "\n",
    "        l: np.array, Lower triangle result (nxn)\n",
    "        u: np.array, Upper triangle result (nxn)\n",
    "        p: np.array, Permutation: 交换后的行顺序 (n)\n",
    "            p = None if sequence=True\n",
    "\n",
    "    Raises:\n",
    "        Exception: 存在为零的主元素\n",
    "    \"\"\"\n",
    "    a = np.array(a, dtype=np.float)  # copy\n",
    "\n",
    "    assert a.shape[0] == a.shape[1]\n",
    "    n = a.shape[0]\n",
    "\n",
    "    if not sequence:\n",
    "        # p 记录行交换的过程，使用「列主元高斯消元法」才使用，否则为 None\n",
    "        p = np.array([k for k in range(n)])\n",
    "        # swap_times:  行交换次数\n",
    "        if isinstance(swap_times, dict):\n",
    "            swap_times['swap_times'] = 0\n",
    "    else:\n",
    "        p = None\n",
    "\n",
    "    for k in range(n-1):\n",
    "        if not sequence:\n",
    "            i_max = k + np.argmax(np.abs(a[k:n, k]))\n",
    "\n",
    "            if i_max != k:\n",
    "                a[[i_max, k]] = a[[k, i_max]]  # swap rows\n",
    "                p[[i_max, k]] = p[[k, i_max]]  # record\n",
    "                swap_times['swap_times'] += 1\n",
    "\n",
    "        if a[k][k] == 0:\n",
    "            raise Exception(\"存在为零的主元素\")\n",
    "\n",
    "        for i in range(k+1, n):\n",
    "            a[i][k] /= a[k][k]  # L @ 严格下三角\n",
    "            for j in range(k+1, n):\n",
    "                a[i][j] -= a[i][k] * a[k][j]  # U @ 上三角\n",
    "\n",
    "    return np.tril(a, k=-1) + np.identity(a.shape[0]), np.triu(a), p\n",
    "\n",
    "\n",
    "def solve_lu(b, l, u, p=None):\n",
    "    \"\"\"用 lu(a) 得到的 `pa=lu` 分解的结果求解原方程组 `ax=b` 的解 x。\n",
    "\n",
    "    若 p 不为 None 则使用「列主元高斯消元」，p 为 None表示使用「顺序高斯消元」。\n",
    "\n",
    "        # `@` means matrix multiplication, refer: https://docs.python.org/reference/expressions.html#binary-arithmetic-operations\n",
    "        b = p @ b if p != None\n",
    "        l @ y = b\n",
    "        u @ x = y\n",
    "\n",
    "    Args:\n",
    "        b: np_array_like, 原方程组的右端常数（n）\n",
    "        l: np_array_like, Lower triangle of lu_seq(a)\n",
    "        u: np_array_like, Upper triangle of lu_seq(a)\n",
    "        p: np_array_like, LU分解中交换后的行顺序\n",
    "            default p=None: 未做行交换，即使用顺序高斯消去法\n",
    "\n",
    "        使用列主元高斯消元法时，l, u, p 使用 lu(a) 得到的结果即可：\n",
    "            solve_lu(b, *lu(a))\n",
    "        或者使用顺序高斯消元：\n",
    "            solve_lu(b, *lu(a, sequence=True))  # p=None\n",
    "\n",
    "    Returns:\n",
    "        x : np.array `ax=b` 的解（n）\n",
    "    \"\"\"\n",
    "    assert np.shape(l) == np.shape(u)\n",
    "    assert np.shape(l)[0] == np.shape(b)[0]\n",
    "\n",
    "    n = np.shape(l)[0]\n",
    "\n",
    "    # do swap\n",
    "    if p is not None:\n",
    "        b = [b[v] for v in p]\n",
    "\n",
    "    # L * y = b\n",
    "    y = np.zeros(n, dtype=np.float)\n",
    "    y[0] = b[0]\n",
    "    for i in range(1, n):\n",
    "        bi = b[i]\n",
    "        for j in range(0, i):\n",
    "            bi -= y[j] * l[i][j]\n",
    "        y[i] = bi / l[i][i]\n",
    "    # print(y)\n",
    "\n",
    "    # U * x = y\n",
    "    x = np.zeros(n, dtype=np.float)\n",
    "    x[n-1] = y[n-1] / u[n-1][n-1]\n",
    "    for i in range(n-2, -1, -1):  # from n-2 (included) to 0 (included)\n",
    "        yi = y[i]\n",
    "        for j in range(i+1, n):\n",
    "            yi -= x[j] * u[i][j]\n",
    "        x[i] = yi / u[i][i]\n",
    "    # print(x)\n",
    "\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse_iteration_lu(A, m0=1, u0=None, eps=1e-8, max_steps=500, verbose=False):\n",
    "    \"\"\"反幂法（inverse iteration a.k.a. inverse power method）计算矩阵特征值、特征向量\n",
    "    \n",
    "    基于 LU 分解\n",
    "    \n",
    "    Args:\n",
    "        A:   np_array_like 待求特征值的矩阵 (nxn)\n",
    "        m0:  float 初始特征值\n",
    "                default m0=1\n",
    "        u0: np_array_like 初始特征向量（n）：要求无穷范式=1，通常取 u0 = (1, 1, ..., 1)\n",
    "                default u0=None: 取 u0 = (1, 1, ..., 1)\n",
    "        eps: float 精度要求\n",
    "                default eps=1e-8\n",
    "        max_steps: int 最大迭代次数\n",
    "                default max_steps=1000\n",
    "        verbose: bool, 若为 True 则打印出每一步的结果\n",
    "                default verbose=False\n",
    "        \n",
    "    Returns:\n",
    "        (m, u, k): 在 max_steps 次迭代以内得到第一组的满足 eps 的结果\n",
    "            m: float 所求主特征值\n",
    "            u: np.array 相应的特征向量\n",
    "            k: int 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 u0 尺寸不匹配\n",
    "                    或 u0 = 0\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "        \n",
    "    _shape = np.shape(A)\n",
    "    if len(_shape) < 2 or _shape[0] != _shape[1]:\n",
    "        raise ValueError(f\"unexpected A, shape: {_shape}\")\n",
    "    \n",
    "    if u0 is not None:\n",
    "        if len(u0) != _shape[0]:\n",
    "            raise ValueError(f\"A (shape={_shape}) and u0 (len={len(u0)}) not match\")\n",
    "        if np.all(u0 == 0):\n",
    "            raise ValueError(f'bad u0: u0 == 0')\n",
    "    else:  # not u0\n",
    "        u0 = np.ones(_shape[0])\n",
    "        \n",
    "    m = m0\n",
    "    u = u0\n",
    "    \n",
    "    lupA = lu(A)\n",
    "\n",
    "    for k in range(int(max_steps)):\n",
    "        if verbose:\n",
    "            print(k, 1/m, u)\n",
    "        \n",
    "        m_prev = m\n",
    "        \n",
    "        v = solve_lu(u, *lupA)\n",
    "        mi = np.argmax(np.abs(v))\n",
    "        m = v[mi]\n",
    "        u = v / m\n",
    "                \n",
    "        if abs(m - m_prev) <= eps:\n",
    "            break\n",
    "        \n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). \"\n",
    "                        f\"The last result: 1/m = {1/m}, u={u}\")\n",
    "    \n",
    "    if verbose:\n",
    "        print('result of inverse_iteration_lu:', 1/m, u, k+1)\n",
    "        \n",
    "    return 1/m, u, k+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1.0 [1. 1. 1.]\n",
      "1 2.7499999999999996 [-0.25        0.41666667  1.        ]\n",
      "2 -1.7368421052631577 [-0.39473684 -0.25877193  1.        ]\n",
      "3 -1.5223065250379362 [-0.27587253 -0.28244815  1.        ]\n",
      "4 -1.773137541073161 [-0.24854257 -0.33791831  1.        ]\n",
      "5 -1.8507138028093015 [-0.22932389 -0.36045702  1.        ]\n",
      "6 -1.906734835399633 [-0.21874548 -0.37508306  1.        ]\n",
      "7 -1.939421265284303 [-0.21209945 -0.38385382  1.        ]\n",
      "8 -1.9604633754301677 [-0.20791023 -0.38945545  1.        ]\n",
      "9 -1.9739760929204886 [-0.20520426 -0.39306055  1.        ]\n",
      "10 -1.9828015632828078 [-0.20343978 -0.3954137   1.        ]\n",
      "11 -1.9885994355225236 [-0.2022801  -0.39695986  1.        ]\n",
      "12 -1.992428450842125 [-0.20151431 -0.39798092  1.        ]\n",
      "13 -1.994964998382758 [-0.201007   -0.39865733  1.        ]\n",
      "14 -1.9966489581948166 [-0.20067021 -0.39910639  1.        ]\n",
      "15 -1.997768464464223 [-0.20044631 -0.39940492  1.        ]\n",
      "16 -1.998513415490028 [-0.20029732 -0.39960358  1.        ]\n",
      "17 -1.9990094345024974 [-0.20019811 -0.39973585  1.        ]\n",
      "18 -1.9993398409805128 [-0.20013203 -0.39982396  1.        ]\n",
      "19 -1.9995599908120016 [-0.200088   -0.39988266  1.        ]\n",
      "20 -1.9997067035591083 [-0.20005866 -0.39992179  1.        ]\n",
      "21 -1.9998044881537032 [-0.2000391  -0.39994786  1.        ]\n",
      "22 -1.999869667263002 [-0.20002607 -0.39996524  1.        ]\n",
      "23 -1.9999131152833085 [-0.20001738 -0.39997683  1.        ]\n",
      "result of inverse_iteration_lu: -1.9999420785330355 [-0.20001158 -0.39998455  1.        ] 24\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-1.9999420785330355, array([-0.20001158, -0.39998455,  1.        ]), 24)"
      ]
     },
     "execution_count": 256,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "inverse_iteration_lu(A, eps=1e-5, max_steps=40, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "测试一下二者的运行速度："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 测试运行速度\n",
    "\n",
    "import time\n",
    "\n",
    "def time_cost(func):\n",
    "    \"\"\"print time cost of a function func\"\"\"\n",
    "    def wrapper(*args, **kwargs):\n",
    "        s = time.time()\n",
    "        ret = func(*args, **kwargs)\n",
    "        e = time.time() - s\n",
    "        print(f'{func.__name__} executed in {e: .4f} s')\n",
    "        return ret\n",
    "    return wrapper\n",
    "\n",
    "eps = 1e-5\n",
    "max_steps = 1e6\n",
    "\n",
    "@time_cost\n",
    "def inverse_iter_np_solve(A):\n",
    "    return inverse_iteration(A, eps=eps, max_steps=max_steps)\n",
    "\n",
    "@time_cost\n",
    "def inverse_iter_lu_solve(A):\n",
    "    return inverse_iteration_lu(A, eps=eps, max_steps=max_steps)\n",
    "\n",
    "@time_cost\n",
    "def call_numpy_linalg_eig(A):\n",
    "    return np.linalg.eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "call_numpy_linalg_eig executed in  0.0003 s\n",
      "(-0.034991233156967005+0j)\n",
      "inverse_iter_np_solve executed in  0.0005 s\n",
      "-0.034991232623709595\n",
      "inverse_iter_lu_solve executed in  0.0031 s\n",
      "-0.03499123262370956\n"
     ]
    }
   ],
   "source": [
    "# A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "\n",
    "# 这里是随机的，不一定能算，有时候要多跑几次才能出来正确的结果\n",
    "A = np.random.random((10, 10))\n",
    "while np.linalg.det(A) == 0:\n",
    "    print('na')\n",
    "    A = np.random.random((10, 10))\n",
    "# print(A)\n",
    "\n",
    "r_eig = call_numpy_linalg_eig(A)\n",
    "reals = [x for x in r_eig[0] if not np.iscomplex(x)]\n",
    "i = np.argmin(np.abs(reals))\n",
    "print(reals[i])\n",
    "\n",
    "r_snp = inverse_iter_np_solve(A)\n",
    "print(r_snp[0])\n",
    "\n",
    "r_slu = inverse_iter_lu_solve(A)\n",
    "print(r_slu[0])\n",
    "\n",
    "assert str(r_snp[0])[:6] == str(r_slu[0])[:6], f'r_snp={r_snp[0]}, r_slu={r_slu[0]}'\n",
    "assert str(r_slu[0])[:5] in str(r_eig[0]), f'r_slu ({r_slu[0]}) not in r_eig: {r_eig[0]}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 374,
   "metadata": {},
   "outputs": [],
   "source": [
    "import traceback\n",
    "\n",
    "def origin_translation(A, p, x0=None, eps=1e-8, max_steps=1000, verbose=False):\n",
    "    \"\"\"原点位移算法\n",
    "    (搜不到这个方法英文是什么)\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 待求特征值的矩阵 (nxn)\n",
    "        p:  float 位移因子\n",
    "        x0: np_array_like, 初始向量（n）, 要求 max(abs(x0)) == 1\n",
    "                default x0=None: 取 x0 = (1, 1, ..., 1)\n",
    "        eps: float, 控制精度\n",
    "                default eps=1e-8\n",
    "        max_steps: 最大迭代次数\n",
    "                default max_steps=1000\n",
    "        verbose: 打印出每步的结果，default verbose=False\n",
    "    \n",
    "    Returns:\n",
    "        (eig, x, k)\n",
    "            eig: float, A 的靠近 p 的特征值\n",
    "            x: np.array, 单位化的特征向量\n",
    "            k: int, 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 x0 尺寸不匹配\n",
    "                    或 x0 = 0\n",
    "                    或 max(u0) != 1\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "        \n",
    "    _shape = np.shape(A)\n",
    "    if len(_shape) < 2 or _shape[0] != _shape[1]:\n",
    "        raise ValueError(f\"unexpected A, shape: {_shape}\")\n",
    "    \n",
    "    if x0 is not None:\n",
    "        if len(x0) != _shape[0]:\n",
    "            raise ValueError(f\"A (shape={_shape}) and x0 (len={len(x0)}) not match\")\n",
    "        if np.all(x0 == 0):\n",
    "            raise ValueError(f'bad x0: x0 == 0')\n",
    "    else:  # not u0\n",
    "        x0 = np.ones(_shape[0])\n",
    "        \n",
    "    _I = np.identity(_shape[0])\n",
    "    \n",
    "    eig = p\n",
    "    x = x0\n",
    "    \n",
    "    for k in range(int(max_steps)):\n",
    "        if verbose:\n",
    "            print(k, eig, x)\n",
    "        \n",
    "        eig_prev = eig\n",
    "        \n",
    "        try:\n",
    "            x = np.linalg.solve(A - eig * _I, x)\n",
    "        except np.linalg.LinAlgError as e:\n",
    "            # traceback.print_exc()\n",
    "            x = np.linalg.solve(A - eig * _I + eps, x)\n",
    "                \n",
    "        q = np.argmax(np.abs(x))\n",
    "        xq = x[q]\n",
    "        \n",
    "        eig = eig + 1 / xq\n",
    "\n",
    "        \n",
    "        x = x / xq\n",
    "        \n",
    "        if np.all(np.abs(eig - eig_prev) <= eps):\n",
    "            break\n",
    "    else:\n",
    "        raise Exception(f\"cannot reach eps ({eps}) after max_steps ({max_steps}). \"\n",
    "                                f\"The last result: eig = {eig}, x={x}\")\n",
    "        \n",
    "    if verbose:\n",
    "        print('result of origin_translation:', eig, x, k+1)\n",
    "        \n",
    "    return eig, x, k+1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 316,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 10 [1. 1. 1.]\n",
      "1 10.709090909090909 [0.47272727 1.         0.74545455]\n",
      "2 11.005485349227762 [0.50058964 1.         0.74989723]\n",
      "3 11.000001815769398 [0.50000023 1.         0.74999988]\n",
      "4 11.000000000000188 [0.5  1.   0.75]\n",
      "result of origin_translation: 11.0 [0.5  1.   0.75] 5\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(11.0, array([0.5 , 1.  , 0.75]), 5)"
      ]
     },
     "execution_count": 316,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "origin_translation(A, 10, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 302,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerating_max_val_eig(A, x0=None, eps_pi=1e-1, eps_ot=1e-8, max_steps_pi=1000, max_steps_ot=1000, verbose=False):\n",
    "    \"\"\"动态原点位移算法求按模最大特征值\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 待求特征值的矩阵 (nxn)\n",
    "        x0: np_array_like, 初始向量（n）, 要求 max(abs(x0)) == 1\n",
    "                default x0=None: 取 x0 = (1, 1, ..., 1)\n",
    "        eps_pi: float,「正幂法」迭代的控制精度\n",
    "                default eps_pi=1e-1\n",
    "        eps_ot: float,「原点位移」算法的控制精度，必须满足 eps > eps_ot\n",
    "                default eps_ot=1e-8\n",
    "        max_steps_pi: int,「正幂法」迭代的最大迭代次数\n",
    "                default max_steps_pi=1000\n",
    "        max_steps_ot: int,「原点位移」算法的最大迭代次数\n",
    "                default max_steps_ot=1000\n",
    "        verbose: bool:\n",
    "                - False: 不打印任何过程信息;\n",
    "                - True: 打印出「正幂法」、「原点位移」的完整迭代过程;\n",
    "                default verbose=False\n",
    "    \n",
    "    Returns:\n",
    "        (eig, eigv, k)\n",
    "            eig: float, A 的按模最大特征值\n",
    "            eigv: np.array, 单位化的特征向量\n",
    "            k: 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 x0 尺寸不匹配\n",
    "                    或 x0 = 0\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "    # 异常输入检测会在 power_iteration 中完成，此处不必重复辛劳\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"    power_iteration:\")\n",
    "    \n",
    "    eig, eigv, k_pi = power_iteration(A, u0=x0, eps=eps_pi, max_steps=max_steps_pi, verbose=verbose)\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"    origin_translation:\")\n",
    "    \n",
    "    eig, eigv, k_ot = origin_translation(A, eig, x0=eigv, eps=eps_ot, max_steps=max_steps_ot, verbose=verbose)\n",
    "    \n",
    "    return eig, eigv, k_pi + k_ot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 303,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    power_iteration:\n",
      "0 1 [1. 1. 1.]\n",
      "1 17.0 [0.41176471 1.         0.58823529]\n",
      "2 9.470588235294118 [0.52795031 1.         0.82608696]\n",
      "3 11.58385093167702 [0.49276139 1.         0.72600536]\n",
      "4 10.831635388739945 [0.50200485 1.         0.75743775]\n",
      "5 11.049799514875502 [0.49945569 1.         0.74783731]\n",
      "result of power_iteration: 10.985906091381928 [0.50014864 1.         0.75061668] 6\n",
      "    origin_translation:\n",
      "0 10.985906091381928 [0.50014864 1.         0.75061668]\n",
      "1 10.999996033673042 [0.49999985 1.         0.74999938]\n",
      "2 11.000000000001119 [0.5  1.   0.75]\n",
      "result of origin_translation: 11.0 [0.5  1.   0.75] 3\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(11.0, array([0.5 , 1.  , 0.75]), 9)"
      ]
     },
     "execution_count": 303,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "accelerating_max_val_eig(A, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 377,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerating_min_val_eig(A, x0=None, eps_ii=1e-1, eps_ot=1e-8, max_steps_ii=1000, max_steps_ot=1000, verbose=False):\n",
    "    \"\"\"动态原点位移算法求按模最小特征值\n",
    "    \n",
    "    Args:\n",
    "        A:  np_array_like, 待求特征值的矩阵 (nxn)\n",
    "        x0: np_array_like, 初始向量（n）, 要求 max(abs(x0)) == 1\n",
    "                default x0=None: 取 x0 = (1, 1, ..., 1)\n",
    "        eps_ii: float,「反幂法」迭代的控制精度\n",
    "                default eps_ii=1e-1\n",
    "        eps_ot: float,「原点位移」算法的控制精度，必须满足 eps > eps_ot\n",
    "                default eps_ot=1e-8\n",
    "        max_steps_ii: int,「反幂法」迭代的最大迭代次数\n",
    "                default max_steps_ii=1000\n",
    "        max_steps_ot: int,「原点位移」算法的最大迭代次数\n",
    "                default max_steps_ot=1000\n",
    "        verbose: bool:\n",
    "                - False: 不打印任何过程信息;\n",
    "                - True: 打印出「正幂法」、「原点位移」的完整迭代过程;\n",
    "                default verbose=False\n",
    "    \n",
    "    Returns:\n",
    "        (eig, eigv, k)\n",
    "            eig: float, A 的按模最大特征值\n",
    "            eigv: np.array, 单位化的特征向量\n",
    "            k: 迭代次数\n",
    "        \n",
    "    Raises:\n",
    "        ValueError: 参数 A 不是方阵，\n",
    "                    或 A 和给定 x0 尺寸不匹配\n",
    "                    或 x0 = 0\n",
    "        Exception:  无法在max_steps 次迭代以内得到满足精度 eps 的结果\n",
    "    \"\"\"\n",
    "    # 异常输入检测会在 power_iteration 中完成，此处不必重复辛劳\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"    inverse_iteration:\")\n",
    "    \n",
    "    eig, eigv, k_ii = inverse_iteration(A, u0=x0, eps=eps_ii, max_steps=max_steps_ii, verbose=verbose)\n",
    "    \n",
    "    if verbose:\n",
    "        print(\"    origin_translation:\")\n",
    "    \n",
    "    eig, eigv, k_ot = origin_translation(A, eig, x0=eigv, eps=eps_ot, max_steps=max_steps_ot, verbose=verbose)\n",
    "    \n",
    "    return eig, eigv, k_ii + k_ot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    inverse_iteration:\n",
      "0 1.0 [1. 1. 1.]\n",
      "1 2.75 [-0.25        0.41666667  1.        ]\n",
      "2 -1.7368421052631582 [-0.39473684 -0.25877193  1.        ]\n",
      "result of inverse_iteration: -1.522306525037936 [-0.27587253 -0.28244815  1.        ] 3\n",
      "    origin_translation:\n",
      "0 -1.522306525037936 [-0.27587253 -0.28244815  1.        ]\n",
      "1 -1.8991931095340484 [-0.22075327 -0.37282221  1.        ]\n",
      "2 -1.9913582399414025 [-0.20172412 -0.39769765  1.        ]\n",
      "3 -1.9999267327390964 [-0.20001466 -0.39998046  1.        ]\n",
      "4 -1.9999999946318703 [-0.2 -0.4  1. ]\n",
      "result of origin_translation: -2.0 [-0.2 -0.4  1. ] 5\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(-2.0, array([-0.2, -0.4,  1. ]), 8)"
      ]
     },
     "execution_count": 378,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "accelerating_min_val_eig(A, verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "寻找一个较优的 eps_{pi, ii}:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 306,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "max:\n",
      "0 8\n",
      "1 9\n",
      "2 11\n",
      "3 12\n",
      "4 14\n",
      "5 15\n",
      "min:\n",
      "0 7\n",
      "1 8\n",
      "2 11\n",
      "3 15\n",
      "4 21\n",
      "5 26\n"
     ]
    }
   ],
   "source": [
    "A = [[2, 3, 2], [10, 3, 4], [3, 6, 1]]\n",
    "\n",
    "print(\"max:\")\n",
    "for i in range(6):\n",
    "    eps_i = 10 ** (-i)\n",
    "    _, _, k = accelerating_max_val_eig(A, eps_pi=eps_i)\n",
    "    print(i, k)\n",
    "    \n",
    "print(\"min:\")\n",
    "for i in range(6):\n",
    "    eps_i = 10 ** (-i)\n",
    "    _, _, k = accelerating_min_val_eig(A, eps_ii=eps_i)\n",
    "    print(i, k)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对这个问题 eps_{pi, ii} 取 1 就比较好。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验内容题目"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 376,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "-- 取初值 u0 = [1, 1, 1]\n",
      "\n",
      "正幂法：(5.999999284744433, array([ 1.        , -0.99999982,  0.99999982]), 24)\n",
      "加速法：(6.000000333333467, array([ 1.        , -0.99999947,  0.99999987]), 20)\n",
      "\n",
      "反幂法：(1.0000004180226194, array([4.18175751e-07, 1.00000000e+00, 9.99999582e-01]), 13)\n",
      "加速法(1.0000000000000009, array([-1.00994543e-15,  1.00000000e+00,  1.00000000e+00]), 9)\n",
      "\n",
      "-- 取初值 u0 = [2.001, 1.999, 0]\n",
      "\n",
      "正幂法：(5.999999476043831, array([ 1.        , -0.99999987,  0.99999987]), 35)\n",
      "加速法：(5.999999666666533, array([ 1.        , -0.99999947,  0.99999987]), 30)\n",
      "\n",
      "反幂法：(1.0000004184290925, array([4.18429144e-07, 1.00000000e+00, 9.99999582e-01]), 14)\n",
      "加速法(1.0000000000000009, array([-9.81794114e-16,  1.00000000e+00,  1.00000000e+00]), 10)\n",
      "\n",
      "-- 取初值 u0 = [0, 1, 1]\n",
      "\n",
      "正幂法：(1, array([0., 1., 1.]), 1)\n",
      "加速法：(0.9999999999999999, array([0., 1., 1.]), 4)\n",
      "\n",
      "反幂法：(1.0, array([2.77555756e-17, 1.00000000e+00, 1.00000000e+00]), 1)\n",
      "加速法(0.9999999999999999, array([0., 1., 1.]), 4)\n",
      "\n",
      "-- 2.5 附近特征值：\n",
      "原点位移： (3.0, array([ 1. ,  0.5, -0.5]), 5)\n"
     ]
    }
   ],
   "source": [
    "A = [[4, -1, 1],\n",
    "     [-1, 3, -2],\n",
    "     [1, -2, 3]]\n",
    "\n",
    "for u0 in [[1, 1, 1], [2.001, 1.999, 0], [0, 1, 1]]:\n",
    "    print(f\"\\n-- 取初值 u0 = {u0}\")\n",
    "\n",
    "    print(\"\\n正幂法：\", end='')\n",
    "    maxlp = power_iteration(A, u0=u0, eps=1e-6)\n",
    "    print(maxlp)\n",
    "\n",
    "    print(\"加速法：\", end='')\n",
    "    maxla = accelerating_max_val_eig(A, x0=u0, eps_pi=1e-3, eps_ot=1e-6)\n",
    "    print(maxla)\n",
    "\n",
    "    print(\"\\n反幂法：\", end='')\n",
    "    minli = inverse_iteration(A, u0=u0, eps=1e-6)\n",
    "    print(minli)\n",
    "\n",
    "    print(\"加速法\", end='')\n",
    "    minla = accelerating_min_val_eig(A, x0=u0, eps_ii=1e-3, eps_ot=1e-6)\n",
    "    print(minla)\n",
    "\n",
    "print(\"\\n-- 2.5 附近特征值：\")\n",
    "r = origin_translation(A, 2.5)\n",
    "print(\"原点位移：\", r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
